Optimized Gaussian exponents for Goedecker-Teter-Hutter pseudopotentials 
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We have optimized the exponents of Gaussian i and p basis functions for the elements H, B-F, and Al-Cl 
using the pseudopotentials of Goedecker, Teter, and Hutter [Phys. Rev. B 54, 1703 (1996)] by minimizing 
the total energy of dimers. We found that this procedure causes the Gaussian to be somewhat more localized 
than the usual procedure, where the exponents are optimized for atoms. We further found that three exponents, 
equal for s and p orbitals, are sufficient to reasonably describe the electronic structure of all elements that we 
have studied. For Li and Be results are presented for pseudopotentials of Hartwigsen et al. [Phys. Rev. B 58, 
3641 (1998)]. We expect that our exponents will be useful for density functional theory studies where speed is 
important. 

PACS numbers: 71.15.Ap,31.15.E- 



I. INTRODUCTION 

Density functional theorji^ (DFT) has become a standard 
tool for electronic structure calculations in condensed-matter 
physics. Using computer clusters, nowadays even solids with 
several hundreds of atoms per unit cell can routinely be stud- 
ied. For molecular dynamics simulations of large systems, 
however, existing codes are still too slow to be able to per- 
form large numbers of time steps within a reasonable amount 
of time. To overcome this problem, faster DFT codes are 
needed. To this purpose two approximations are very useful. 
First of all, pseudopotentials eliminate the need to take into 
account core electrons and simplify the computation for the 
valence electrons by being softer than the original all-electron 
potentials, thus reducing the number of basis functions that 
have to be included in practical calculations. Secondly, the 
use of localized (atom-centered) basis states makes it possi- 
ble to compute many quantities, such as, the electronic den- 
sity, by means of procedures that scale linearly with the num- 
ber of atoms in the unit cell [so-called order(A^) methods]. 
Unfortunately, although many different pseudopotentials have 
been published throughout the years, the problem of finding 
the best possible local basis functions for these pseudopoten- 
tials is only rarely addressed. One reason is that the form 
of such optimized basis states depends on the pseudopoten- 
tials used and is thus not universal. Another reason is perhaps 
that local basis states, especially Gaussians, are mainly used 
by quantum chemists, who are often interested in getting the 
most accurate results possible, which arguably speaks against 
approximating the all-electron potential by pseudopotentials. 
We would here like to stress that our goal is to perform cal- 
culations on large systems, trading accuracy for speed, and to 
compete with tight-binding methods predominantly used for 
such systemsi^i^ 

Toward our proposed goal we present, in this paper, 
optimized Gaussian exponents for the pseudopotentials of 
Goedecker, Teter, and Hutteri^ We chose Gaussians as basis 
functions in part because a lot of literature exists that de- 
scribes how to calculate Hamiltonian matrix elements in an 
efficient, analytical way, and in part because the tails of Gaus- 
sian functions decay very rapidly. The latter property can be 



exploited to speed up the ensuing computer code. We have se- 
lected the pseudopotentials of Goedecker, Teter, and Hutter^^ 
for the following reasons: (1) They have an analytical form, 
which is characterized by only a few parameters. (2) Their 
form has been optimized for efficient real-space DFT com- 
putations. (3) They are available for the first two rows of 
the periodic table. (4) pseudopotentials of a slightly gener- 
alized form are available for elements throughout the peri- 
odic table [See Hartwigsen et al. (Ref. 7)]. (5) Their validity 
has been thoroughly tested,^ among other ways by comparing 
bond lengths of 34 small molecules with all-electron DFT re- 
sults. (6) They can be used in ABINIT,** a freely available DFT 
program, which uses a plane wave basis. We used ABINIT to 
obtain highly accurate reference results for dimers of the ele- 
ments that we studied as well as for a number of other small 
molecules. The remainder of this paper is organized as fol- 
lows. In Sec, nil we describe our DFT calculations. Our op- 
timized exponents are presented in Sec.|IIIl We conclude in 
SecHVl 



II. COMPUTATIONAL PROCEDURE 

For our DFT calculations we wrote a program following the 
work of Obara and Saika,'^ which gives recursive expressions 
for overlap, kinetic, and Hartree matrix elements of s- and 
p-type orbitals, starting from the much easier expressions for 
.s-type orbitals. Corresponding expressions for the Goedecker- 
Teter-Hutter pseudopotentials were straightforwardly derived. 
The exchange and correlation potential (we used the local den- 
sity approximation of Perdew and Wang'") was calculated on 
a fine mesh of grid points, without further approximations. 
For the Hartree matrix elements the exact four-center integrals 
were calculated. No intermediate fitting procedure was used 
for the charge density. In the selfconsistent loop the density 
matrix, describing the electronic charge density, was mixed 
using Broyden mixing. ' ' Electronic occupation numbers were 
calculated according to the Fermi-Dirac statistics with an elec- 
tronic temperature of 1 mHa. Because our main interest is in 
optimizing exponents for speed, we limited ourselves to Gaus- 
sian basis states with s and p character only, an approximation 
which should work best for light elements. 



As mentioned above, reference results were obtained with 
the computer program ABINIT*^ This program uses periodic 
boundary conditions and a plane wave basis. For the ABINIT 
results presented below, a careful convergence study was 
made with respect to the box size and the number of plane 
waves included to ensure that the total free energies obtained 
were converged to within 1 mHa or better 



III. RESULTS 

Usually Gaussian exponents are optimized for atoms, 
which may not give optimal results for molecules or solids. 
To improve on this situation we used dimers for our optimiza- 
tions. We found that the Gaussian functions became some- 
what more localized as compared to the atomic optimizations, 
which is favorable if speed is important (fewer matrix ele- 
ments need to be calculated). For Si we additionally opti- 
mized the exponents for three atoms arranged on the corners 
of a triangle and for four atoms located on the vertices of a 
square. For these latter optimizations we found no systematic 
change in the exponents. Therefore, we decided to optimize 
our Gaussian basis functions using dimers. 

In our optimization procedure the total free energy of the 
dimers was minimized including different numbers of s and p 
orbitals in the basis. The following configurations were tested: 
Isp, Islp, ls2p, 2sp, 2s\p, 2s2p, 2s3p, 3sp, 3slp, 3s2p, 3s3p, 
Asp, 4s3p, AsAp, where nsmp means that n i-type and m p- 
type Gaussians were used, and nsp indicates that every one 
of n exponents were used for s and p Gaussians. The possi- 
bility of equal s- and /^-exponents was included, because it is 
computationally advantageous (see, for example, Ref. 9). Af- 
ter each optimization, we computed the binding energy of the 
dimer for a wide range of atomic separations. 

Surprisingly, we found that the 3sp data gave close to op- 
timal results for the binding energies of all dimers. The re- 
sulting curves are shown in Figs.[T]and|2] For hydrogen, the 
Asp curve (not shown) was slightly better than the 3sp curve, 
but the energy differences between both curves were small 
compared to the errors that we made in other elements by re- 
stricting ourselves to a basis containing only s and p orbitals. 
For lithium and beryllium the pseudopotentials of Goedecker, 
Teter, and Hutter^ treat the \s electrons as valence electrons, 
and offer therefore only a limited advantage in speed as com- 
pared to all-electron calculations. For this reason, for these 
two elements we decided to use the above-mentioned pseu- 
dopotentials of Hartwigsen et al.,'' which for Li and Be happen 
to be parameterized by the same analytical form as the pseu- 
dopotentials of Goedecker, Teter, and Hutter, so that no ad- 
ditional programming was necessary for these elements. The 
Li and Be pseudopotentials of Hartwigsen et al— that we have 
used treat the \s electrons as core electrons, and should thus 
provide a considerable speed-up as compared to all-electron 
calculations. 

As one can see in Fig. [T]our results for H, Li, Be, B, and C 
are excellent. At first sight, the agreement of our dimer curve 
with the corresponding ABINIT reference curve looks slightly 
worse for Be2 [Fig.fljc)] than for the other dimers [Figs.fTfa), 
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Element 


rf(ao) 






Exponents (aq 






H 


L45 


4, 


.8336700 


0.7962986 


0, 


,1703991 


Li 


5.099 


0, 


,7107769 


0.1151854 


0, 


,0381934 


Be 


4.515 


1, 


,3420164 


0.2543896 


0, 


,0760345 


B 


3.00 


2, 


,2292786 


0.4922699 


0, 


,1321358 


C 


2.35 


3, 


,3582735 


0.7697578 


0, 


,2000381 


N 


2.07 


4, 


,6545372 


1.0441223 


0, 


,2624989 


O 


2.27 


6, 


,3165398 


1.4495363 


0, 


,3579517 


F 


2.61 


8, 


,3059235 


1.9605670 


0, 


,4596180 


Al 


4.64 


0, 


,9496758 


0.2017059 


0, 


,0659053 


Si 


4.29 


1, 


,1985394 


0.2972473 


0, 


,0941132 


P 


3.57 


1, 


,4859314 


0.4104534 


0, 


,1277851 


S 


3.57 


1, 


,7834663 


0.5430840 


0, 


,1651222 


CI 


3.74 


2, 


,1074139 


0.7324463 


0, 


,2130882 



TABLE I: Optimized Gaussian basis set exponents for H, Li-F, and 
Al-Cl. The same exponents should be used for Gaussian s and p 
orbitals. The second column gives the atomic separations id) of the 
dimers for which the exponents were optimized. It should be stressed 
that the exponents for Li and Be were optimized using pseudopoten- 
tials of Hartwigsen et al.J whereas all other exponents were opti- 
mized for Goedecker-Teter-Hutter pseudopotentials.— 



Molecule 


Bond Length (ao) 


Binding 


Energy (Ha) 




isp data 


ABINIT 


'isp data 


ABINIT 


BH 


2.383 


2.366 


-0.1860 


-0.1857 


FH 


1.784 


1.763 


-0.2970 


-0.3026 


AlH 


3.210 


3.161 


-0.1573 


-0.1599 


HCl 


2.451 


2.435 


-0.2212 


-0.2279 


CO 


2.162 


2.127 


-0.5322 


-0.5715 


CS 


3.003 


2.884 


-0.3422 


-0.3912 


PN 


2.930 


2.789 


-0.3978 


-0.4528 


SiO 


2.977 


2.833 


-0.3854 


-0.4368 



TABLE II: Bond lengths and binding energies of diatomic molecules 
obtained using the 3.sp Gaussian basis sets from Table|l] Reference 
results were calculated using the computer program ABINIT. 



[nb),[Tfd), andlHe)]. The reason is that the binding energy of 
Be2 at its equilibrium distance is approximately an order of 
magnitude less than that of, for example, H2 or C2. As a re- 
sult, the 3sp basis set errors look large for Be2 [Fig. [He)], even 
though the absolute error is less than 9 mHa for all distances 
shown (2.5 mHa at the minimum), which should be acceptable 
for most applications. In Fig.[T] starting from nitrogen, the in- 
completeness of our basis set becomes gradually noticed, and 
as a result the agreement between our results and the ABINIT 
reference curves becomes gradually worse. The same behav- 
ior can be noticed in the dimer curves for the second-row el- 
ements, which are shown in Fig.|2] Whereas the 3sp basis set 
gives very reasonable results for Al, Si, and P, the curves for 
sulphur and chlorine differ increasingly from the ABINIT ref- 
erence curves. It would, of course, depend on the application, 
whether the deviations are acceptable. For future reference, 
all our optimized 3sp exponents for H, Li-F, and Al-Cl are 
presented in Table U 

To test, how well the exponents of Table U perform in sys- 
tems other than the dimers for which they were optimized, 
we calculated the binding energies and the bond lengths for 
8 diatomic molecules. Our results are given in Table HIl The 
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FIG. 1 : (Color online) The binding energy as a function of atomic distance for the dimers (a) H2 and (b)-(h) Li2-F2. The curves labeled ABINIT 
show reference results obtained with the computer program ABINIT. The curves labeled "S^p data" show the results that we have obtained 
using optimized basis sets with 3 different exponents, but identical for s and p Gaussians. It should be stressed that the results for (b) Li and 
(c) Be were obtained using pseudopotentials of Hartwigsen et al.r whereas all other curves were calculated using Goedecker-Teter-Hutter 
pseudopotentials.— 
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FIG. 2: (Color online) The binding energy as a function of atomic distance for dimers of the second-row elements Al-Cl. The curves labeled 
ABINIT show reference results obtained with the computer program ABINIT. The curves labeled "3sp data" show the results that we have 
obtained using optimized basis sets with 3 different exponents, but identical for i- and p Gaussians. 



first thing to be noticed is that the bond lengths obtained with 
the Gaussian 3sp basis sets are always larger than the ABINIT 
reference bond lengths, which are also given in Table HIl and 
that the 3sp binding energies are always smaller than the refer- 
ence results, with the exception of BH, for which the 3sp and 
the reference binding energies agree almost perfectly. These 
trends are consistent with what can be observed for the dimers 
(Figs. [Hand |2]l. A more quantitative analysis shows that for 
BH, FH, AlH, HCl, and CO the Gaussian 3sp basis reproduces 
the ABINIT bond lengths to within 2%, and that for CS, PN, 
and SiO the deviations of the bond lengths are of the order 
of 5%. The errors in the binding energies are between 0.3 
and 6.7 mHa for BH, FH, AlH, and CIH, and between 39.3 



and 55 mHa for CO, CS, PN, and SiO, but never more than 
12.5% of the ABINIT reference binding energy. All these er- 
rors lie within the range of what one might reasonably have 
expected from the dimer curves (Figs. [T] and |2]i, and we con- 
clude that our Gaussian 3sp basis sets perform equally well for 
the dimers and the small molecules that we have studied. We 
once again remind the reader that our main goal in the present 
work is to compete with tight-binding theory rather than to 
produce very accurate results for a variety of small molecules 
composed of first and second row elements. 
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IV. CONCLUSION 

In this paper we have optimized Gaussian exponents for 
the Goedecker-Teter-Hutterpseudopotentials'' for all elements 
in the first two rows of the periodic table, with the excep- 
tion of the noble gases and Na and Mg. We decided not 
to include the noble gases in our study, because the local 
density approximation'" is known to perform quite badly for 
these elements. Sodium and magnesium were kept out of 
the scope of this study, because the Goedecker-Teter-Hutter 
pseudopotentials^ treat the 2s and 2p electrons as valence elec- 
trons, which would make practical applications unnecessarily 
expensive in terms of computer time. 

The main conclusion of the present work is that for all ele- 
ments studied it is sufficient to use a Gaussian basis containing 



only s and p orbitals, and furthermore that only three Gaussian 
exponents, which can be chosen equal for the s and the p chan- 
nels, are needed to obtain a reasonable accuracy, as has been 
discussed in the body of this paper. 

We believe that our results constitute an important contribu- 
tion towards the development of fast DFT computer programs 
using both pseudopotentials and Gaussian basis sets. 
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